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Abstract 


Reduction  of  computational  error  is  a  key  issue  in  computing  Lagrangian  trajectories 
using  gridded  velocities.  Computational  accuracy  enhances  from  using  the  first  tenn  (constant 
velocity  scheme),  the  first  two  terms  (linear  uncoupled  scheme),  the  first  three  terms  (linear 
coupled  scheme),  to  using  all  the  four  terms  (nonlinear  coupled  scheme)  of  the  two-dimensional 
interpolation.  A  unified  ‘analytical  form’  is  presented  in  this  study  for  different  truncations. 
Ordinary  differential  equations  for  predicting  Lagrangian  trajectory  are  linear  using  the  constant 
velocity/linear  uncoupled  schemes  (both  commonly  used  in  atmospheric  and  ocean  modeling), 
linear  coupled  scheme  and  nonlinear  using  the  nonlinear  coupled  scheme  (both  proposed  in  this 
paper).  Location  of  the  Lagrangian  drifter  inside  the  grid  cell  is  determined  by  two  algebraic 
equations,  which  are  solved  explicitly  with  the  constant  velocity/linear  uncoupled  schemes,  and 
implicitly  using  the  Newton-Raphson  iteration  method  with  the  linear/nonlinear  coupled 
schemes.  The  analytical  Stommel  ocean  model  on  the /-plane  is  used  to  illustrate  great  accuracy 
improvement  from  keeping  the  first-term  to  keeping  all  the  terms  of  the  two-dimensional 
interpolation. 
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1.  Introduction 


Oceanic  and  atmospheric  motion  can  be  represented  by  Eulerian  and  Lagrangian 
viewpoints.  The  fonner  gives  time-dependent  three-dimensional  (Eulerian)  fields  of  velocity, 
temperature,  salinity,  and  other  variables,  which  are  commonly  represented  in  satellite 
observations,  modeling,  simulation,  and  prediction  at  numerical  grid  points.  The  latter  provides 
continually  changing  characteristics  (temperature,  salinity,  velocity,  etc.)  along  the  fluid 
particles’  trajectories  (i.e.,  Lagrangian  trajectories),  which  are  commonly  represented  in  in-situ 
oceanographic  measurements  by  Argo  floats,  drifters,  and  gliders.  Employing  the  Lagrangian 
trajectories,  water  masses  can  also  be  distinguished  in  terms  of  origin  and/or  destination  and  be 
traced  (Vries  and  Doos  2001).  The  two  types  of  velocity  are  convertible.  Routine  ocean  data 
assimilation  systems  (Galanis  et  al.  2006;  Lozano  et  al.  1996;  Song  and  Colberg  2011;  and  Sun 
1999)  and  data  analysis  methods  such  as  optimal  interpolation  (OI)  (Gandin  1965)  and  optimal 
spectral  decomposition  (OSD)  (Chu  et  al.  2003  a,  b),  can  be  used  for  converting  Lagrangian 
drifter  data  into  gridded  Eulerian-type  data,  and  evaluating  ocean  models  (e.g.,  Chu  et  al.  2001, 
2004).  Several  new  phenomena  were  discovered  after  the  conversion.  Lor  example,  with  the 
OSD  method  new  signals  have  been  identified  such  as  fall-winter  recurrence  of  current  reversal 
from  westward  to  eastward  on  the  Texas-Louisiana  continental  shelf  from  near-surface  drifting 
buoy  and  current-meter  (Chu  et  al.  2005),  and  propagation  of  long  baroclinic  Rossby  waves  at 
mid-depth  (around  1,000  m  deep)  in  the  tropical  north  Atlantic  from  the  Argo  floats  (Chu  et  al. 
2007). 

Consider  water  particles  flowing  with  ocean  currents  in  three-dimensional  space  (x,  y,  z) 
and  time  t,  discretized  into  grid  cells  with  the  spacing  of  (Ax,  Ay,  Az)  and  time  step  of  At,  with 
the  discrete  Eulerian  velocity  filed  represented  by 
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v(X .  y,  ,zk,t,)  =  [u(xt ,  yj  ,zk,t, ),  v(x. ,  y.  ,zk,t,),  w(xi ,  yy ,  z, ,  t,) 


(la) 


Here,  the  subscripts  (/,  j,  k,  /)  represent  the  spatial  and  temporal  discretization.  The  superscsript 
‘A’  means  the  Eulerian  gridded  fields.  Common  interpolation  methods  can  be  used  to  get  four 
dimensional  continuous  velocity  field  from  gridded  field  (la), 

v(x, y,  z,  t )  =  [u(x, y,  z,  t ),  v(x, y,  z,  t ),  w(x, y,z,t )] .  (lb) 

The  position  of  each  fluid  particle,  R(t)  =  [x(f),  y(t),  z(/)],  is  specified  in  the  Lagrangian  system. 
The  connection  between  the  Eulerian  and  Lagrangian  approaches  leads  to  the  ordinary 
differential  equations, 


dx(t ) 
dt 


u(x,  y,  z,  t). 


dy(t ) 
dt 


v(x,y,z,t), 


dz{t ) 
dt 


w(x,y,z,t), 


(2a) 


which  determines  the  trajectory  of  the  particle  if  the  position  is  specified  at  some  initial  instant  in 
its  path  history.  Such  calculation  has  also  been  used  as  the  semi-Lagrangian  scheme  in  ocean 
numerical  modeling  (e.g.,  Chu  and  Fan  2010).  Thus,  the  interpolation  (lb)  is  the  key  in 
calculating  Lagrangian  trajectories  from  gridded  velocity  fields.  For  steady  gridded  velocity 
fields,  the  analytical  solution  exists  for  the  Lagrangian  trajectory  (2a)  inside  one  grid  cell  with 
(lb)  a  highly  truncated  linear  interpolation  in  space  (see  Section  2  for  explanation)  (Doos  1995; 
Blanke  and  Raynaud  1997).  Follow-up  research  has  been  extended  from  steady  to  unsteady 
velocity  fields  with  the  Lagrangian  trajectories  being  calculated  from  time-varying  gridded 
velocity  fields  (Vries  and  Doos  2001). 

Two  sources  of  uncertainty  exist  in  determining  the  Lagrangian  trajectories  from  Eulerian 
flow  field:  (a)  the  knowledge  of  the  smoothness;  and  (b)  the  error  in  the  integration  of  the 
ordinary  differential  equations  (2a).  There  is  a  need  to  estimate  uncertainties  due  to  the  limited 
knowledge  of  the  Eulerian  velocity  (see  Section  2).  To  illustrate  this  point,  consider  the  case  that 
we  only  have  access  to  the  average  of  the  velocity  in  a  cell.  In  this  case  the  trajectories  within  the 
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cell  are  straight  lines,  called  the  constant  velocity  (CV)  scheme.  With  more  knowledge  about  the 
Eulerian  velocity,  for  example,  Vries  and  Doos  (2001)  used  low  order  truncation  in  spatial 
interpolation  [see  (5a)  and  (5b)  in  Section  2]  to  simplify  [u(x,  y,  z,  t),  v(x,  y,  z,  ?)]  in  (2a)  by 
dx(t )  dy(t ) 

— )—  =  Lx(x,t)  =  a0  +  aJ  +  (ax+aJ)x,  — =  L2(y,t)  =  /30  +  /?T  +  (A  +  fi2t)y,  (2b) 
dt  dt 

where  u  depends  on  (x,  t)  only  and  v  depends  on  (y,  t)  only.  Such  a  treatment  leads  to  the 
existence  of  analytical  solutions.  The  following  coefficients  in  (2b)  vanish 

&2  —  Cl^  —  P2  ~  (d^  —  0, 

when  the  Eulerian  flow  field  is  steady.  The  two  functions  in  (2b)  are  represented  by 

A  (x,  t )  =  Zj  (x)  =  a0+  axx ,  L2  ( y ,  t)  =  L2  (y)  =  j30+  fdxy,  (2c) 

In  reality,  for  a  2D  Eulerian  flow  field,  the  velocity  components  u{x,  y,  t),  and  v(x,  y,  t)  in  (2a) 
are  not  necessarily  taken  the  functions  (L\,  L2)  given  by  (2b).  Questions  arise:  What  is  the 
Lagrangian  trajectory  if  the  Eulerian  velocity  components  ( u ,  v)  depend  on  (x,  y)  [more 
realistic]?  Is  there  any  improvement  with  such  a  change?  In  other  words,  what  is  the 
improvement  for  a  steady  Eulerian  flow  field  if  L,(x)  is  changed  into  u(x,  y)  and  L2(y)  is 
changed  into  v(x,  y)?  What  is  the  improvement  for  an  unsteady  Eulerian  flow  field  if  L\(x,  t)  is 
changed  into  u(x,  y,  t)  and  L2(y,  t )  is  changed  into  v(x,  y,  t)  for  an  unsteady  Eulerian  flow  field? 
To  show  the  accuracy  progressive  in  calculation  of  Lagrangian  trajectories,  a  systematical 
analysis  is  presented  in  this  study  for  a  steady  Eulerian  flow  field,  and  will  be  presented  in 
another  paper  in  the  near  future  for  an  unsteady  Eulerian  flow  field.  Division  of  steady  and 
unsteady  Eularian  flow  fields  is  due  to  the  mathematical  complexity. 

The  rest  of  the  paper  is  outlined  as  follows.  Section  2  describes  the  establishment  of 
continuous  velocity  inside  a  grid  cell.  Section  3  depicts  the  calculation  of  Lagrangian  trajectory 
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from  side  to  side  of  a  grid  cell.  Section  4  shows  the  identification  of  starting  grid  cell.  Section  5 
describes  the  Lagrangian  trajectory  across  the  grid  cell.  Section  6  introduces  the  Stommel  ocean 
model  for  the  evaluation.  Section  7  shows  the  accuracy  progressive  from  high  to  no  truncation  of 
the  two-dimensional  interpolation.  Section  8  presents  the  conclusions. 

2.  Establishment  of  Continuous  Velocity  Inside  a  Gridded  Cell 
For  simplicity  without  loss  of  generality,  a  steady-state  two-dimensional  gridded  data  is 
considered.  Let  the  water  particle  be  located  at  (called  a  starting  point,  not  necessary  at  the  grid 
point)  R0  =  (x0,y0)  inside  the  grid  cell:  [x(  <  x0  <  x/+1  y;  <  y0  <  yJ+l  ]  and  let  it  move  using  the 

gridded  data.  Due  to  spatial  variability  of  the  gridded  velocity  data,  the  Lagrangian  velocity 
changes  with  time  although  the  Eulerian  flow  is  steady.  Let  the  velocity  be  given  at  the  four 
comer  points  of  the  grid  cell,  Ftj,  Fi+\,j,  F^+i,  Fi+\j+\.  Here,  F  represents  (u,  v).  For  a  two- 
dimensional  interpolation,  the  velocities  inside  the  ij  grid  cell  can  be  given  by  the  corner  points, 

E  (x,  y )  =  a0  +  a,  (x  -  xM  )  +  a2(y-  y,_, )  +  a3  (x  -  xM  )  (y  -  yH  )  .  (3) 

Let  the  Lagrangian  drifter  travel  from  (xo,  yo)  to  (xi,  y\)  with  the  travel  time  of  r  (Fig.  1),  and  the 
Lagrangian  velocity  components  w(x,  y)  and  v(x,  y)  be  represented  by 


u(x,y)  u(x0y0)  -  ^  ^  ^  ^(x  x0), 

Xj-Xq 

(4a) 

,  ^  x  v{xl,yl)-v(x0,y0) 

v(x,y)  v(x0,y0)=  (y  y0j- 

(4b) 

Substitution  of  (3)  into  (4a)  and  (4b)  leads  to 
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u(x,y)-u(x0y0) 

«(Wo) 


<{x^yj+x) 


Sx 


l_xl-xQ 


1- 


ZlZZo 

Sy 


+  “(*i+l>Jo)- 


x,  —  X, 


r 


j 


l- 


Ji"Jo 


A 


Sx 


Ji"Jo 


Sy 


+  U 


(*<+i>Jy+i) 


Sx  \  Sy  j 

*i“*o  Ji~Jo 


Sx  Sy 


—  k(*o>Jo) 


m(*/+i.Jo)-“(*q.Jo)  ,  M(:yo»J7+i)~“(:yo»Jo)  Ji~Jo  | 


Ax 


Sy 


Ji“Jo 


AxAy 


vAxy0 


r 

V  Ax  Jo 


“  (*0  >  Jo  )  -  : U  [XM  5  Jo  )  -  W  (*0>  yj+l  )  +  U  (*,+l >  J j+1  ) 


(x-x0) 


f— l 

v  A J  y  o  Ai  —  Ao 


Ji-Jo  +  Ji-Jo 


Xj  -  x0  AxAy 


Ax  Ay  _  Ax  Sy 
Ax  Ay 


Ax  Ay 


*i+l  J 


-u 


i,j+ 1 


Ax  Ay 
Ax  Ay 


+ 


/+1J+1 


Ax  Sy 
Ax  Ay 


(x-x0) 


Ji-Jo 


r  Au^ 

v  Ay  Jo  *i-*0 


^  .2  A 

An 


v(x,j)-v(x0,j0)  = 


^Ava 

v  Ay  y0 


+ 


v  AxAy 
^  Av^ 


(Ji-Jo) 


Jo 


Xx-XQ 


(x-x0), 


v Axy0  y,-yQ 


+ 


'  A2v 

.^AyJo 


(xi-x0) 


(j-Jo)  = 


where 


f AF) 

.f(x/+i,Jo)-F(x0,Jo) 

f  AF  ^ 

V  Ax  J 

o 

1 

81 

l  Ay  J 

F(w;+i)-F(w0) 


f  a 2f  a 

V  AxA_y  jQ 


FiJ-FMJ  -FiJ+l+Fi+l,j+l 

Ax  Ay 


Sy 


,  Ax  =  xi+1  -x0,  Sy  =  yj+l-y0, 


(x-x0) 


(5a) 


(5b) 


are  given  from  the  gridded  velocities  as  well  as  the  starting  velocity  [n(xo,  yo),  v(xo,  Vo)]  with  the 
starting  position  (xo,  yo)-  Vries  and  Doos  (2001)  only  keep  the  first  term  in  the  bracket  of  the 
right  hand  side  of  each  equation  in  (5a)  and  (5b)  and  argued  that  inclusion  of  last  two  terms  was 
impossible  to  give  a  general  analytical  solution  although  it  may  be  important  in  the  case  of 
strongly  curved  streamlines. 
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Eqs.(5a)  and  (5b)  can  be  rewritten  into  a  more  general  form, 
u(x,y)  =  u(xo,y0)  +  Ax(x-x0),  v(x,y)  =  v(xo,y0)  + Av(y  -  y0) ,  (6) 


where 


A  = 


A  = 


f  An^ 
\Axj0 

v4 yjo 


yi~y0 


< An > 

vA yj0  Xi-x0 
f  AW 


y{-y0 


\—} 

K^yJo 

f  A2v  ^ 

vAxA y 


(yi-y0)> 


(Xj-Xq). 


7o 


Substitution  of  (6)  into  (2a)  leads  to 
dx(t ) 


dt 


dy(t) 


dt 


(x, y)  =  u(x0, y0 )  +  Ax  (x - x0) , 


=  v(x, y)  =  v(xo,y0)  +  Ay  (y -y0). 


which  have  the  following  solutions, 

x(t)  =  x0+u(x0,y0)t,  y(t)  =  y0+v(x0,y0)t,  ifAx=Ay=0, 


x(t)  =  x0+"{X°;^  (eA*‘-l),  if  4*0, 
(An/ Ax)o 

y(0  =  y0  +  (eA/  - 1),  if  Av*  0. 


(AW  A y)0 


(7) 


(8a) 

(8b) 


(9) 

(10a) 


(10b) 


The  solutions  (9),  (10a),  and  (10b)  imply  that  the  Lagrangian  drifter  never  moves  if  the  starting 
velocity  equals  zero,  i.e.,  no  =  n(x o,  >’o)  =  0,  and  Vo  =  v(xo,  Vo)  =  0. 

For  a  sufficiently  short  travel  time  t  with  the  Lagrangian  drifter  still  being  inside  the  ij 
grid  cell,  the  location  (xi,  y\)  can  be  easily  obtained  if  Ax  =0,  Ay  =  0,  or  (Ax  ,  Ay)  are  given  [i.e., 
keeping  the  first  term  in  the  right-hand  side  of  (7)], 

Xj  =x0+n(x0,y0)r,  y1  =y0  +v(x0,y0)r,  if  Ax  =0,Ay  =0,  (11a) 
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x1=x0  + 


u(x0,y0) 
(A m  /  Ax)o 


04r-l),  V  =Ao  + 


v(Wo)  ,  v 


(Av/A y\ 


(ey  ~1),  if  4  = 


V  Ax  y0 


.4-  = 


^Ava 
vAf  70 


■  (Hb) 


For  more  general  cases  [keeping  the  first  two  or  all  terms  in  the  right-hand  side  of  (7)],  the 
location  (xi,  y\)  satisfies  the  following  two  non-linear  algebraic  equations, 

«(Wo) 


Xj  =  x0  + 


(Aw  /  Ax)( 


-  {exp  [4(x15  jA] r  - 1}  > 


(12a) 


yi=y°+l^i\y)  (12b) 

which  are  solved  by  the  Newton-Raphson  iteration  method. 

3.  Lagrangian  Trajectory  from  Side  to  Side  of  a  Grid  Cell 

Various  truncation  of  (3)  leads  to  accuracy  increase  in  calculating  the  Lagrangian 
trajectory  (inside  the  ij  grid  cell)  from  the  gridded  velocities  at  the  four  corners  of  the  ij  grid  cell. 
If  only  the  first  tenn  in  the  right-hand  side  of  (3)  is  used,  i.e.,  Ax  =  0,  Ay  =  0,  the  two  ordinary 
differential  equations  (8a)  and  (8b)  become 


whose  solutions  are 


dx(t )  dy(t) 

— j~A  =  u(x0,y0),  — -^  =  v(x0,y0): 
cit  cit 


x(t)  =  x o  +u(x0,y0)t,  y(t)=y0+v(x0,y0)t, 


(13) 


(14) 


which  is  called  the  constant  velocity  (CV)  scheme  since  the  velocity  components  [m(xo,  Vo), 
v(xo,  yo)]  are  constant  during  the  movement  of  the  Lagarangian  drifter  inside  the  ij  grid  cell. 

If  the  first  two  terms  in  the  right-hand  side  of  (3)  are  used,  i.e., 

^  Alt  ^ 


A  = 


V  Ax  y0 


A  = 


^Ava 

vAfyo 


(15) 
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the  two  differential  equations  (8a)  and  (8b)  do  not  depend  on  (xi,  v i )  and  have  analytical 
solutions  (Doos  1995;  Blanke  and  Raynaud  1997;  Vries  and  Doos  2001), 


x{t)  =  x0  + 


“(x0,y0) 

(Am/  Ax)q 


(e  -l),  y(t)  =  y0  + 


y(w.) 

(Av/Ay)0 


At 


(e  y  - 1)  . 


(16) 


It  is  called  the  linear  uncoupled  (LUC)  method. 

If  the  first  three  terms  in  the  right-hand  side  of  (3)  are  used,  i.e., 


4(Wi): 


'  Am2 
vAxJ0 


Li  ~y0 


rAu^ 

\Ayj0  xi~xo 


Y(Wi): 


^Av^  YAv"! 

V  Ax  J 


vAyy0 


oLi-Lo 


(17) 


the  two  differential  equations  (8a)  and  (8b)  depend  on  (xi,  y\),  which  represents  the  end  point  of 
the  trajectory.  It  is  called  the  linear  coupled  (LC)  scheme  since  the  two  velocity  components 
[m(xi,  v i ),  v(xi,  v i )]  depend  on  both  xi  andyi  linearly.  If  all  the  four  terms  in  the  right-hand  side 
of  (3)  are  used,  i.e., 


Y(Wi)  = 


AJxl,y1)  = 


Am 

V  Ax  y0 


^AvA 
vAf  7o 


Li-Lo 


Au 

\Ayj0  xi~xo 


+ 


Am 

AxA_y 


(u  ~y0), 


J  0 


Xi-X0 


f  A2v  ^ 


\Ax)0yx-y 0 


Ax  A y 


(x,  -x0) , 


(18a) 


(18b) 


J  o 


the  two  differential  equations  (8a)  and  (8b)  also  depend  on  (xi,  v i ) .  It  is  called  the  nonlinear 
coupled  (NLC)  scheme  since  the  two  velocity  components  [n(x i,  y i),  v(xi,  yi)]  depend  on  both  xi 
and  y\  nonlinearly.  During  the  integration  of  (8a)  and  (8b),  the  location  (xi,  yi)  is  determined 
from  solving  the  two  nonlinear  algebraic  equations  (12a)  and  (12b)  using  the  Newton-Raphson 
iteration  method. 


4.  Identification  of  Starting  Grid  Cell 

Let  a  Lagrangian  trajectory  start  from  the  initial  location  (xoo,  >’oo).  If 

x,  <  x00  <  X/+, ,  (19a) 
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yj<yoo<yj+i’ 


(19b) 


the  point  (x0o,  yoo)  is  located  inside  the  ij  grid  cell.  As  the  trajectory  hits  the  side  or  corner  of  the 
initial  grid  cell  at  the  location  (xo,  yo),  it  is  important  to  determine  which  the  next  grid  cell  is  for 
the  advance  of  the  trajectory.  The  location  (xo,  yo)  is  called  the  starting  point  of  the  next  grid  cell. 
The  Lagrangian  trajectory  is  always  calculated  across  the  grid  cell  from  (xo,yo)  at  the  left  or  right 
side  (Fig.  2),  the  upper  or  lower  side  (Fig.  3),  and  the  grid  point  (Fig.  4). 

Let  Fig.  2  be  taken  as  an  example  for  the  illustration  since  Fig.  3  is  similar  but  in  the  y- 
direction.  For  n0  ^  0  (Fig.  2a),  the  point  (xo,  yo)  is  located  at  the  left  (right)  side  and  will  move 
to  the  right  (left)  grid  cell  if  n0  >  0  (no  <  0).  For  n0  =  0  and  v0  0  (Fig.  2b  and  2c), 
determination  of  the  next  grid  cell  depends  on  both  signs  of  [  v0,  (An  /  Ay)0  ].  Solutions  (9),  (10a), 
and  (10b)  require  one  component  of  (no,  v0)  non-zero.  For  no  =  0,  vo  must  be  non-zero.  With  vo  > 
0,  the  starting  point  (xo,  yo)  is  located  in  the  right  cell  for  (An  /  Ay)0  >  0 ,  and  in  the  left  cell  for 

(An  /  Ay)0  <  0  (Fig.  2b).  With  v0  <  0,  the  starting  point  (xo,  yo)  is  located  in  the  right  cell  for 
(An  /  Ay)0  <  0 ,  and  in  the  left  cell  for  (An  /  Ay)0  >  0  (Fig.  2c).  For  n0  =  0  and  v0  =  0,  the 
trajectory  stays  at  (xo,yo)  forever. 

For  (xo,  yo)  located  at  the  corner  of  the  grid  cell  (i.e.,  at  the  grid  point  such  as  at  xo  =  x,-,  yo 
=  Vj  (Fig.  4),  the  point  (xo,  yo)  will  move  to  the  upper  right  cell  for  (no  >0,  v0  >  0),  the  upper 
left  cell  for  (no  <0,  v0  >  0),  the  lower  left  cell  for  (no  <0,  v0  <  0),  and  the  lower  right  cell  for 
(no  >0,  vo  <  0).  With  v0  =  0,  the  point  (xo,  yo)  will  move  to  the  upper  right  cell  for  [no  >  0, 
(A v/  Ax)0  >  0  ],  the  upper  left  cell  for  [no  <  0,  (Av/  Ax)0  <  0],  the  lower  left  cell  for  [no  <  0, 
(A v/ Ax)0  >  0  ],  and  the  lower  right  cell  for  [n0  >  0,  (Av/ Ax)0  <  0].  With  no  =  0,  the  point  (xo, 
yo)  will  move  to  the  upper  right  cell  for  [vo  >  0,  (An  /  Ay)0  >  0  ],  the  upper  left  cell  for  [vo  >  0, 
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(A u  /  A y)0  <  0],  the  lower  left  cell  for  [v0  <  0,  (An  /  Ay)0  >  0  ],  and  the  lower  right  cell  for  [v0  < 
0,  (Au  /  Ay)0  <  0]. 

5.  Lagrangian  Trajectory  Across  Grid  Cell 

The  solutions  (9)  or  (10a,  b)  are  valid  within  a  given  grid  cell.  If  (xi,  y i)  hits  the  corner 
point  or  side  of  the  grid  cell  (x&,  yb),  i.e.,  (xi  =  xb,  yi  =  }’/,),  this  ending  point  (xb,  v/>)  is  treated  as 
the  starting  point  for  the  next  grid  cell,  and  determined  by  the  travel  time  in  the  x-direction 


T,.  =  < 


Xb~X  0 

<(w<>)' 


AAxb^b) 


In 


(x.-Xp^x,,^)  |  1 
u  (xo’To) 


for  Ax  =  0,  Av  =  0 


for  A,_  ^  0,  A,.  ^  0 


(20a) 


andy-direction, 


T„  =  < 


yb-y0 

v{x^y*Y 


In 


Ay(Xb^yb) 


(yh-y0)Av(xh,yh)  1 1 
v(x0>y0) 


for  A  =  0,  A  =  0 


for  A,_^0,A  ^0 


(20b) 


The  Lagrangian  trajectory  hits  the  comer  of  the  ij  grid  cell  [i.e.,  one  of  the  four  grid  points  (x„ 
yj),  (xi+ 1,  vy),  (xj+i,  yj),  (x,+i,  y,+i)]  ifrx  =  r  =  r  .  The  Lagrangian  trajectory  hits  the  side  of  the  cell 

if  rx  ^  ry .  For  tx  >  Ty ,  it  hits  the  upper  side  if  Vo  >  0  and  hits  the  lower  side  if  v0  <  0.  For 
tx  <  ry ,  it  hits  the  right  side  if  uo  >  0  and  hits  the  left  side  if  no  <  0  (Fig.  5). 

For  the  Lagrangian  trajectory  hitting  the  cell  side,  one  of  Xb  and  y/,  takes  the  grid  point 
location  (one  of  x,-,  x,-+i,  yj,  yr\  )  and  the  other  is  obtained  from  solving  an  algebraic  equation  with 
the  constraint  of  rx  =  r  , 


v(xo,y0)(xb -xo)-u{x0’y0)(yb -y0)  =  0’ 


(21a) 
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for  the  CV  scheme, 


(xb-x0)  Ax  (x0,y0)  (yb-y0)  Ay  (x0 ,  y0 ) 


'(Wo) 


;(Wo) 


=  exp 


4r(Wo) 

^(Wo) 


(21b) 


for  the  LUC  scheme, 


Ay(xb’yb)1 


n 


(xb-x0)Ax(xb,yh)  |  { 
u{x0,y0) 


AAxb^yb) 1 


n 


(yh-y0)Ay(xb^yb) 

v{x0,y0) 


=0,  (21c) 


for  the  LC  and  NLC  schemes.  Since  one  of  (xb,  yb)  is  given,  (21a)  and  (21b)  are  linear  algebraic 
equations,  which  is  solved  easily  and  explicitly.  However,  (21c)  is  a  nonlinear  algebraic 
equation,  which  is  solved  using  the  Newton-Raphson  method.  After  (x^  >’b)  being  obtained,  the 
travel  time  r  is  detennined,  and  in  turn  the  Lagragian  trajectory  is  obtained  before  hitting  the 
grid  cell  side  (or  corner)  using  (10a)  or  (10b)  for  0  <t  <t  (i.e.,  dashed  curve  in  Fig.  1). 

It  is  also  noted  that  during  the  integration  the  velocity  components  ( u ,  v)  are  set  to  zero 
under  the  conditions, 


u  =  0  if 

u 

<  10  ^s1. 

hi 

o 

ll 

> 

V 

Ar 

Ay 

<10 


-10-1 


(22) 


The  relative  displacement  components  (|x/  Ax|,  | y  I  Ay|  )  are  rounded  with  the  accuracy  of  10 


-9 


6.  Stommel  Ocean  Model  on  the  /-Plane 

Stommel  (1948)  designed  an  ocean  model  to  explain  the  westward  intensification  of  wind 
driven  ocean  currents.  Consider  a  rectangular  ocean  with  the  origin  of  a  Cartesian  coordinate 
system  at  the  southwest  corner  (Fig.  6).  The  x*-  and  y*-  axes  point  eastward  and  northward, 
respectively.  Here,  the  superscript  denotes  dimensional  variables.  The  boundaries  of  the 
ocean  are  at  x*  =  0,  L;  and  y*  =  0,  b.  The  ocean  is  considered  as  a  homogeneous  and 
incompressible  layer  of  constant  depth  D  when  at  rest.  When  currents  occur  as  in  the  real  ocean, 
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the  depth  differs  from  D  everywhere  by  a  small  perturbation.  Due  to  the  incompressibility,  a 
streamfunction  i//*  is  defined  by 


dxF* 


u*  =  — 


dy: 


v*  = 


fir 


where  u*  and  v*  are  components  of  the  velocity  vector  in  the  x*  and  y*  directions.  The  surface 
wind  stress  is  taken  as  -F  cos(ny/b).  The  component  frictional  forces  are  taken  as  -Ru  and  -Rv, 
where  R  is  the  frictional  coefficient.  The  Coriolis  parameter /is  also  introduced.  For  a  constant/ 
an  equation  was  derived  for  the  streamfunction  y/*, 


f  d2 


32  A 


fic  *2  dy  *2 


*  •  <ny  \ 

Y*  -  ~Y  sm( — - — )  . 

b 


(23) 


where  y  =  Fn  /  (Rb) .  The  rigid  boundary  conditions  are  given  by 

Y(  0,  y*)  =  y*)  =  0)  =  b)  =  0. 

The  independent  and  dependent  variables  are  non-dimensionalized  by 
x  =  x*/A-0.5,  y  =  y*/ b  —  0.5,  y/  =  y/*K2  / (yb2) . 


(24) 


(25) 


For  simplicity  without  loss  of  generality,  the  dimensional  parameters  (2,  b)  are  chosen  such  as 
b  =  l.  The  analytical  solution  of  Eq.(23)  in  the  non-dimensional  form  is  given  by  (Fig.  6b) 


y/{x,y)  =  sm(ny) 


'  l-e-1  x 

1 - 

e-e 


e-\ 


e-e 


(26) 


with  0<x<l,  0  <  y  <  \  and  the  maximum  value, 


w  =1-2 

t  max 


eV2-e-m 


max  -I 

e  —  e 


(27) 


The  non-dimensional  velocity  components  of  the  Stommel  model  (its,  vs)  are  given  by 
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us(x,y)  =  =  - 7rcos(7ry ) 

vs(x,y)  =  ^  =  -sm(ny) 
ox 


1- 


\-e 


-e  —■ 


e  —  e 


l-e-1 


e-l 


-e  -- 


e-e 


e  —  e 


e-l 
e  —  e1 

A 


y 


(28) 


7.  Accuracy  Progressive  among  the  Four  Schemes 


The  non-dimensional  ocean  basin  is  discretized  by  Ax  =  Ay  =  0.02  .  The  velocity  components 
are  calculated  at  the  grid  points  ( u,j ,  v,j)  (i  =  1,  2,  51;  j  =1,  2,  51)  using  (28).  With  the 

gridded  Eulerian  velocity  fields  {uup  Vy),  the  continuous  velocity  fields  [w(x,  y),  v(x,  y)]  are 
obtained  using  (6)  with  four  different  methods  (CV,  LUC,  LC,  and  NLC).  Since  the  Stommel 
model  on  the  f-plane  is  symmetric  (Fig.  6)  the  initial  location  is  selected  by 


*00  =0.14,  Too  =0-0-  (29) 

which  is  2.5  times  away  from  the  boundary  than  from  the  center  of  the  circulation.  Eq.(26) 
shows  that  the  stream-function  at  (xo,  yo)  is  given  by 


y/Q  =^(0.14,  0.0)  =  0.1045.  (30) 

Since  the  Stommel  ocean  model  has  the  steady-state  analytical  solution,  the  Lagrangian  drifter  is 
supposed  to  move  along  any  closed  streamline  (Fig.  7a),  which  means  that  the  Lagrangian 
trajectory  coincides  with  the  streamline  and  should  be  closed.  The  two  differential  equations  (8a) 
and  (8b)  are  integrated  using  the  four  schemes  (CV,  LUC,  LC,  NLC)  for  computing  Ax(x,  y)  and 
Ay(x,  v)  with  the  Lagragian  trajectory  moving  around  the  ocean  basin  up  to  100  circles. 

First,  the  analytical  streamline  (Fig.  7a)  is  used  to  evaluate  the  accuracies  of  the  CV, 
LUC,  LC,  and  NLC  schemes  (Figs.  7b-e).  The  Lagrangian  trajectory  is  not  a  closed  circle  using 
the  CV,  LUC,  and  LC  schemes.  Using  the  CV  scheme,  it  hits  the  ocean  boundary  after  8.375 
circles  (Fig.  7b),  using  the  LUC  scheme,  it  hits  the  ocean  boundary  after  26  circles  (Fig.  7c), 
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using  the  LC  scheme,  it  does  not  hit  the  ocean  boundary  after  100  circles  (but  not  a  closed 
streamline  with  v|/-value  changing  to  0.032  after  100  circles)  (Fig.  7d).  However,  using  the  NLC 
scheme,  it  is  exactly  the  same  as  the  analytical  streamline  after  100  circles  (\|/ -value  kept  as 
0.1045)  (Fig.  7e).  Since  vp  =  0  at  the  lateral  boundary,  the  following  criterion 

H<KT6,  (31) 

is  used  to  identify  the  Lagragian  trajectory  hitting  the  lateral  boundary. 

Second,  the  initial  streamfunction  <//0  is  used  to  evaluate  the  four  methods.  The  smaller 

the  difference  of  the  numerical  v|/-value  against  y/0,  the  more  accurate  the  scheme  would  be.  Fig. 

8  shows  the  dependence  of  the  v|/-value  versus  the  circle  of  the  Lagrangian  trajectory  around  the 
ocean  basin.  The  zero  value  of  the  streamfunction  indicates  the  ocean  boundary.  The  vp -value 
reduces  from  0.1045  to  0  at  the  8.375th  (26th)  circle  using  the  CV  (LUC)  method;  and  to  0.032  at 
the  100th  circle  using  the  LC  method.  The  v|/-value  keeps  0.1045  after  100th  circle  using  the 
NLC  method. 

Third,  the  relative  root  mean  square  error  (RRMSE)  of  the  streamfunction  is  used  to 
evaluate  the  accuracy, 

RRMSE  =  5  (32) 

Wa 

where  y/a  is  the  analytical  streamfunction.  The  comparison  is  conducted  with  three  different 
initial  locations  (Table  1)  with  associated  ^-values  (0.1045,  0.08752,  0.06143)  (Fig.  9).  For  the 
same  initial  location  (jcoo,  yoo),  RRMSE  increases  from  0  to  1.0  in  5-8  circles  using  the  CV 
method,  and  in  17  circles  using  the  LUC  method,  increase  from  0  to  around  0.7  in  100  circles 
using  the  LC  method.  RRMSE  keeps  near  0  in  100  circles  using  the  NLC  scheme  (Fig.  10a).  For 
the  same  method,  RRMSE  increases  as  the  initial  location  changing  towards  the  boundary  (from 
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Case  1  to  Case  3  in  Table  1).  To  further  investigate  the  performance  of  the  NLC  method, 
RRMSE  (in  10'  )  is  plotted  in  one  circle  for  the  three  initial  locations  (Fig.  10b).  The  oscillation 
of  RRMSE  is  noted  with  the  largest  (smallest)  amplitude  for  Case-3  (Case-1).  The  minimum 
RRMSE  values  occur  when  the  trajectory  passing  four  points  (0,  y),  (x-,  0),  (0,  y+),  (x+,  0)  with 
either  u  =  0  or  v  =  0.  The  2-D  calculation  becomes  1-D  calculation,  and  greatly  decreases  the 
RRMSE. 

The  CPU  time  comparison  is  based  either  on  the  first  5  circles  (Table  2)  or  first  100 
circles  (or  hitting  the  boundary)  (Table  3)  of  the  Lagrangian  drifter  around  the  streamline  of  the 
Stommel  model.  Since  the  Stommel  model  is  steady  state,  the  Lagrangian  trajectory  coincides 
with  the  Eulerian  streamline.  The  calculated  Lagrangian  trajectory  has  less  (more)  deviation  to 
the  streamline  using  more  (less)  accurate  scheme  with  accuracy  decreasing  from  the  NLC,  LC, 
LUC,  to  CV  scheme  (see  Fig.  7).  Thus,  the  Lagrangian  drifter  moves  the  shortest  distance  per 
circle  using  the  NLC  scheme  and  longest  distance  using  the  CV  method.  Except  the  CV  scheme 
(consuming  the  least  CPU  time  -  0.0031  s  per  circle),  for  the  first  5  circles,  the  NLC  scheme 
consumes  less  CPU  time  per  circle  (0.03432  s)  than  the  LUC  scheme  (0.06552  s)  and  the  LC 
scheme  (0.03744  s)  (Table  2),  and  up  to  100  circles  (or  hitting  the  boundary),  the  NLC  scheme 
consumes  comparable  CPU  time  per  step  (0.000660  s)  as  the  LUC  scheme  (0.000636  s)  and  the 
LC  scheme  (0.0000646  s).  The  lowest  CPU  time  per  circle  and  per  step  using  the  CV  is  caused 
by  the  simplest  calculation  of  the  Lagrangian  trajectory  [i.e.,  Eq.(9)]. 

8.  Conclusion 

(1)  Two  sources  of  uncertainty  in  determining  the  Lagrangian  trajectories  from  the 
Eulerian  velocity  are  identified:  (a)  the  knowledge  of  the  smoothness;  and  (b)  the  error  in  the 
integration  of  the  ordinary  differential  equations.  This  study  especially  shows  the  process  of 
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establishing  series  of  accuracy  progress  schemes  (CV,  LUC,  LC,  NLC)  with  different  knowledge 
of  smoothness  for  calculating  Lagrangian  trajectory  using  the  gridded  velocity  held  through 
different  truncations  of  a  two-dimensional  interpolation.  Ah  the  four  schemes  are  within  the 
same  analytical  framework  using  two  coefficients  (Ax,  Ay)  with  the  time-dependence  of  the 
Lagrangian  trajectory  analytical:  linear  for  the  CV  scheme,  and  exponential  for  the  rest  schemes 
(LUC,  LC,  NLC). 

(2)  Accuracy  increases  with  the  change  of  the  two  coefficients  (Ax,  Av).  When  Ax  =  Av  =  0, 
the  Lagragian  velocity  components  use  the  starting  velocity  (no,  vo)  (the  CV  scheme),  the 
accuracy  is  the  lowest.  When  (Ax,  A  v)  are  truncated  at  the  first  term  of  the  right-hand  side  in  Eq 
(7),  the  Lagragian  velocity  component  u  depends  on  x,  and  v  depends  on  y  only  (the  LUC 
scheme),  the  accuracy  is  the  lower.  When  (Ax,  A  ,)  are  truncated  at  the  first  two  terms  of  the  right- 
hand  side  in  Eq  (7),  the  Lagragian  velocity  components  (n,  v)  depend  on  (x,  y)  linearly  (the  LC 
scheme),  the  accuracy  is  higher.  When  (Ax,  A ,,)  keep  ah  the  three  terms  of  the  right-hand  side  in 
Eq  (7),  the  Lagragian  velocity  components  ( u ,  v)  depend  on  (x,  y)  nonlinearly  (the  NLC 
scheme),  the  accuracy  is  the  highest.  The  Lagrangian  trajectory  is  obtained  explicitly  using  the 
CV  and  LUC  schemes  and  implicitly  using  the  LC  and  NLC  schemes  with  the  Newton-Raphson 
iteration  method. 

(3)  The  non-dimensional  (length  of  1.0  in  both  x  and  v  directions)  Stommel  ocean  model 
(steady-state  with  analytical  solution)  on  the  /-plane  is  used  for  the  evaluation.  The  Lagrangian 
trajectory  is  calculated  from  the  initial  location  at  the  distance  of  0.14  to  the  center  of  the  ocean 
basin  using  the  four  schemes  (CV,  LUC,  LC,  and  NLC)  from  the  gridded  velocity  data  obtained 
from  the  analytical  Stommel  ocean  model.  The  Lagrangian  trajectory  is  accurately  determined 
with  no  deviation  from  the  steamline  even  after  the  Lagrangian  drifter  moving  around  the  ocean 
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basin  100  circles  using  the  NLC  scheme;  less  accurately  determined  with  deviation  from  the 
streamline  using  the  LC  scheme;  inaccurately  determined  with  evident  deviation  from  the 
streamline  (hitting  the  ocean  boundary  after  26  circles)  using  the  LUC  scheme;  and  very 
inaccurately  determined  with  large  deviation  from  the  streamline  (hitting  the  ocean  boundary 
after  8.375  circles)  using  the  CV  scheme.  The  CV  scheme  consumes  the  least  CPU  time.  The 
NLC  scheme  consumes  comparable  CPU  time  as  the  LUC  and  LC  schemes. 

(4)  High  accuracy  with  no  evident  increase  of  CPU  time  makes  the  NLC  scheme  a 
promising  schemefor  calculating  Lagrangian  trajectory  from  gridded  velocity  data  especially 
with  strongly  curved  streamlines. 

(5)  Calculation  of  Lagrangian  trajectories  from  2D  gridded  velocity  field  described  here 
is  easy  to  extend  to  3D  gridded  velocity  field  by  changing  2D  grid  cell  into  3D  grid  volume.  For 
the  procedure  identified  in  Section  4,  the  trajectory  starts  from  a  surface  (or  grid  point)  of  the 
grid  volume  (xo,  yo,  zo)  rather  than  a  side  (or  grid  point)  of  the  grid  cell  (xo,  yo),  and  ends  at  a 
surface  (or  grid  point)  ( Xb ,  yb,  Zb)  rather  than  a  side  (or  grid  point)  of  the  grid  cell  (x/„  v/i).  The 
starting  point  for  the  next  grid  volume  is  determined  by  equalizing  the  three  travel  times  ( zx ,  zy, 
zz)  [similar  to  (20a),  and  (20b)],  rx  =  ry  =  rz ,  which  provides  two  algebraic  equations  of  (y/,,  z/,) 

for  the  CV  scheme  [similar  to  (21a)],  the  LUC  scheme  [similar  to  (21b)],  and  the  LC  and  NLC 
schemes  [similar  to  (21c)].  The  two  algebraic  equations  are  solved  by  2-D  Newton-Raphson 
method. 

(6)  The  semi-Lagrangian  method  combines  both  Eulerian  and  Lagrangian  points  of  view. 
The  fluid  variable  is  discretized  on  an  Eulerian  grid,  but  is  advanced  in  time  using  the  equation 
similar  to  Eq.(2a).  The  algorithms  in  the  context  of  calculations  of  drifter  trajectories  (e.g.,  the 
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366  CV,  LUC,  LC,  NLC  schemes)  can  be  applied  to  the  semi-Lagrangian  methods  in  ocean 

367  modeling. 

368  (7)  The  limitation  of  this  study  is  that  only  an  analytical  steady  ocean  model,  i.e.,  the 

369  Stommel  model,  is  used  for  evaluating  the  four  schemes.  In  the  context  of  practical  application  to 

370  the  trajectories  of  drifters  driven  by  oceanographic  fields,  it  will  be  useful  to  examine  the 

371  properties  of  these  algorithms  in  the  near  future  under  somewhat  more  realistic  conditions;  for 

372  instance  output  from  an  eddy  resolving  ocean  model  (i.e.,  unsteady  Eulerian  gridded  flow  field). 

373  Acknowledgments: 

374  This  research  was  supported  by  the  Office  of  Naval  Research  (N0001413AF00002)  and  the 

375  Naval  Oceanographic  Office  (N6230612P000123). 

376  References 

377  Blanke,  B.,  and  S.  Raynaud,  1997:  Kinematics  of  the  Pacific  Equatorial  Undercurrent:  An 

378  Eulerian  and  Lagrangian  approach  from  GCM  results.  J.  Phys.  Oceanogr.,  27,  1038-1053. 

379 

380  Chu,  P.C.,  and  C.W.  Fan,  2010:  Space-time  transformation  in  flux-form  semi-Lagrangian 

381  schemes.  Terr.  Atmos.  Ocean.  Sci.,  21  (1),  17-26. 

382 

383  Chu,  P.C.,  S.H.  Lu,  and  Y.C.  Chen,  2001:  Evaluation  of  the  Princeton  Ocean  Model  using  the 

384  South  China  Sea  Monsoon  Experiment  (SCSMEX)  data.  J.  Atmos.  Oceanic  Technol.,  18,  1521- 

385  1539. 

386  Chu,  P.C.,  G.H.  Wang,  and  C.W.  Fan,  2004:  Evaluation  of  the  U.S.  Navy’s  Modular  Ocean  Data 

387  Assimilation  System  (MODAS)  using  the  South  China  Sea  Monsoon  Experiment  (SCSMEX) 

388  data.  J.  Oceanogr.,  60,  1007-1021. 

389  Chu,  P.C.,  L.M.  Ivanov,  T.P.  Korzhova,  T.M.  Margolina,  and  O.M.  Melnichenko,  2003: 

390  Analysis  of  sparse  and  noisy  ocean  current  data  using  flow  decomposition,  Part  1:  Theory.  J. 

391  Atmos.  Oceanic  Technol.,  20,  478-491. 

392  Chu,  P.C.,  L.M.  Ivanov,  T.P.  Korzhova,  T.M.  Margolina,  and  O.M.  Melnichenko,  2003: 

393  Analysis  of  sparse  and  noisy  ocean  current  data  using  flow  decomposition.  Part  2:  Application  to 

394  Eulerian  and  Lagrangian  data.  J.  Atmos.  Oceanic  Technol.,  20,  492-5 12. 

395  Chu,  P.C.,  L.M.  Ivanov,  and  O.M.  Melnichenko,  2005:  Fall-winter  current  reversals  on  the 

396  Texas-Louisiana  continental  shelf,  J.  Phys.  Oceanogr.,  35,902-910. 


20 


397 

398 

399 

400 

401 

402 

403 

404 

405 

406 

407 

408 

409 

410 

411 

412 

413 

414 

415 

416 

417 

418 

419 

420 

421 

422 

423 

424 

425 

426 

427 


Chu,  P.C.,  L.M.  Ivanov,  O.V.  Melnichenko,  and  N.C.  Wells,  2007:  Long  baroclinic  Rossby 
waves  in  the  tropical  North  Atlantic  observed  from  profiling  floats,  J.  Geophvs.  Res.,  112, 
C05032,  doi:  10.1 029/2006JC003698. 

Doos,  K.,  1995:  Interocean  exchange  of  water  masses,  J.  Geophys.  Res.,  100,  13499-13514. 

Gandin,  L.S.,  1965:  Objective  analysis  of  meteorological  fields,  Israel  Program  for  Scientific 
Translation,  Jerusalem,  242  pp. 

Galanis,  G.,  P.  Louka,  P.  Katsafados,  I.  Pytharoulis  and  G.  Kallos,  2006:  Applications  of 
Kalman  filters  based  on  non-linear  functions  to  numerical  weather  predictions,  Ann.  Geophys., 
24,  2451-2460. 

Lozano,  C.J.,  A.R.  Robinson,  H.G.  Arango,  A.  Gangopadhyay,  Q.  Sloan,  P.J.  Haley,  L. 
Anderson,  and  W.  Leslie,  1996:  An  interdisciplinary  ocean  prediction  system:  Assimilation 
strategies  ana  structured  data  models,  Elesevier  Oceanogr.  Series,  61,  413-452. 

Song,  Y.T.,  and  F.  Colberg,  2011:  Deep  ocean  warming  assessed  from  altimeters,  Gravity 
Recovery  and  Climate  Experiment,  in  situ  measurements,  and  a  non-Boussinesq  ocean  general 
circulation  model,  J.  Geophys.  Res.,  116,  C02020,  doi:10.1029/2010JC006601. 

Sun,  L.C.,  1999:  Data  inter-operability  driven  by  oceanic  data  assimilation  needs,  Marine 
Technol.  Soc.  J.,  33  (3),  55-66. 

Stommel,  H.  M.,  1948:  The  westward  intensification  of  wind-driven  ocean  circulation,  Trans. 
Amer.  Geophys.  Union,  29,  202-206. 

Vries,  P.D.,  and  K.  Doos,  2001:  Calculating  Lagrangian  trajectories  using  time-dependent 
velocity  fields,  J.  Atmos.  Oceanic  Technol.,  18,  1092-1101. 


21 


428  Table  1.  Three  initial  locations  and  associated  analytical  \|/u-values. 


429 


Case 

Xo 

To 

430 

431 

1 

0.14 

0.0 

0.1045  432 

433 

2 

0.24 

0.0 

0.08752  434 

435 

3 

0.34 

0.0 

0.06143  436 

437 

438 

439 

440 

441  Table  2.  Comparison  of  CPU  (unit:  s)  for  the  first  5  circles  among  the  four  methods. 


CV 

LUC 

LC 

NLC 

CPU  for  5  circles 

0.0156 

0.3276 

0.  1872 

0.  1716 

CPU  per  circle 

0.0031 

0.06552 

0.03744 

0.03432 

442 

443 

444  Table  3.  Comparison  of  CPU  (unit:  s)  for  the  Lagrangian  trajectories  either  hitting  the  boundary 

445  or  up  to  100  circles  among  the  four  methods. 


CV 

LUC 

LC 

NLC 

Circles  either  hitting 
the  boundary  or  100 

8.375 

26 

100 

100 

Total  CPU 

0.0312 

2.8704 

8.2057 

3.6036 

CPU  per  circle 

0.0037 

0.1104 

0.0821 

0.0360 

Total  steps 

1176 

4517 

12696 

5200 

CPU  per  step 

0.0000265 

0.000636 

0.000646 

0.000660 

446 

447 
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448  Figure  Captions 

449  Fig.  1 .  Illustration  of  a  Lagrangian  trajectory  [x(f),  y(t)\  (dashed  curve)  from  (xo,  yo)  to  (xi,  Vi) 

450  inside  the  ij  grid  cell. 

451  Fig.  2.  Determination  of  initial  grid  cell  with  the  initial  location  of  the  Lagragian  trajectory 

452  (xo,  yo)  located  at  xo  =  x,-  for  (a)  u0  ^  0 ,  (b)  uo  =  0,  v0  >  0,  and  (c)  uq  =  0,  v0  <  0. 

453  Fig.  3.  Determination  of  initial  grid  cell  with  the  initial  location  of  the  Lagragian  trajectory 

454  (xo,  yo)  located  at  yo  =  y,  for  (a)  v0  ^  0 ,  (b)  vo  =  0,  uo  >  0,  and  (c)  vo  =  0,  uo  <  0. 

455  Fig.  4.  Determination  of  initial  grid  cell  with  the  initial  location  of  the  Lagragian  trajectory 

456  (xo,  yo)  located  at  xo  =  x„  yo  =  y,-. 

457  Fig.  5.  Determination  of  the  side  of  the  grid  cell  for  the  Lagrangian  trajectory  crossing. 

458  Fig.  6.  Stommel  ocean  model  on  the/-plane:  (a)  ocean  geometry,  and  (b)  streamfunction  (mVs) 

459  (after  Stommel  1948). 

460  Fig.  7.  Calculated  Lagragian  trajectories  with  the  initial  location  (0.14,  0.00)  and  i|/o  =  0.1045 

461  using  (a)  analytical  solution,  (b)  CV  scheme,  (c)  LUC  scheme,  (d)  LC  scheme,  and  NLC  scheme. 

462  Fig.  8.  Temporal  evolution  of  v|/-values  of  the  Lagrangian  trajectory  calculated  with  the  four 

463  schemes.  It  is  noted  that  the  time  is  represented  by  the  number  of  circles  around  the  ocean  basin. 

464  Fig.  9.  Streamlines  with  three  different  initial  locations:  (0.14,  0.00)  (solid  curve,  vpo  =  0.1045), 

465  (0.24,  0.00)  (dotted  curve,  v|/0  =  0.08752),  and  (0.34,  0.00)  (dashed  curve,  \po  =  0.06143). 

466  Fig.  10.  (a)  Temporal  evolution  of  RRMSE  of  streamfunction  of  the  Lagrangian  trajectory 

467  calculated  with  the  four  schemes  using  the  four  schemes  with  three  different  initial  locations,  (b) 

468  Zoom-in  temporal  evolution  of  RRMSE  of  streamfunction  of  the  Lagrangian  trajectory 

469  calculated  with  the  NLC  scheme  (vertical  scale  is  nearly  three  orders  magnitude  smaller).  It  is 

470  noted  that  the  time  is  represented  by  the  number  of  circles  around  the  ocean  basin. 
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476 

477 

478 

479  Fig.  1 .  Illustration  of  a  Lagrangian  trajectory  [. x{t ),  y{t)]  (dashed  curve)  from  (x0,  yo)  to  (xi,  v  i ) 

480  inside  the  ij  grid  cell. 
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483  Fig.  2.  Determination  of  initial  grid  cell  with  the  initial  location  of  the  Lagragian  trajectory 

484  (xo,  yo)  located  at  xo  =  x,-  for  (a)  u0  ^  0 ,  (b)  uq  =  0,  vq  >  0,  and  (c)  uo  =  0,  vo  <  0. 
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Fig.  3.  Determination  of  initial  grid  cell  with  the  initial  location  of  the  Lagragian  trajectory 
(xo,  yo)  located  at  yo  =  y/  for  (a)  v0  ^  0 ,  (b)  vo  =  0,  uq  >  0,  and  (c)  vo  =  0,  uo  <  0. 
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498  Fig.  4.  Determination  of  initial  grid  cell  with  the  initial  location  of  the  Lagragian  trajectory 

499  (x0,  yo)  located  at  x0  =  xh  y()  =  y7. 
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509  Fig.  5.  Determination  of  the  side  of  the  grid  cell  for  the  Lagrangian  trajectory  crossing. 
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Fig.  6.  Stommel  ocean  model  on  the/-plane:  (a)  ocean  geometry,  and  (b)  streamfimction  (mVs) 
(after  Stommel  1948). 
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518  Fig.  7.  Calculated  Lagragian  trajectories  with  the  initial  location  (0.14,  0.00)  and  i|/o  =  0.1045 

519  using  (a)  analytical  solution,  (b)  CV  scheme,  (c)  LUC  scheme,  (d)  LC  scheme,  and  NLC  scheme. 
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525  Fig.  8.  Temporal  evolution  of  v|/-values  of  the  Lagrangian  trajectory  calculated  with  the  four 

526  schemes.  It  is  noted  that  the  time  is  represented  by  the  number  of  circles  around  the  ocean  basin. 
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531 

532  Fig.  9.  Streamlines  with  three  different  initial  locations:  (0.14,  0.00)  (solid  curve,  v|/0  =  0.1045), 

533  (0.24,  0.00)  (dotted  curve,  vj/0  =  0.08752),  and  (0.34,  0.00)  (dashed  curve,  \|/o  =  0.06143). 
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537  Fig.  10.  (a)  Temporal  evolution  of  RRMSE  of  streamfunction  of  the  Lagrangian  trajectory 

538  calculated  with  the  four  schemes  using  the  four  schemes  with  three  different  initial  locations,  (b) 

539  Zoom-in  temporal  evolution  of  RRMSE  of  streamfunction  of  the  Lagrangian  trajectory 

540  calculated  with  the  NLC  scheme  (vertical  scale  is  nearly  three  orders  magnitude  smaller).  It  is 

541  noted  that  the  time  is  represented  by  the  number  of  circles  around  the  ocean  basin. 
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